Biophysical larval dispersal models of observed bonefish (Albula vulpes) spawning events in Abaco, The Bahamas: An assessment of population connectivity and ocean dynamics

Biophysical models are a powerful tool for assessing population connectivity of marine organisms that broadcast spawn. Albula vulpes is a species of bonefish that is an economically and culturally important sportfish found throughout the Caribbean and that exhibits genetic connectivity among geographically distant populations. We created ontogenetically relevant biophysical models for bonefish larval dispersal based upon multiple observed spawning events in Abaco, The Bahamas in 2013, 2018, and 2019. Biological parameterizations were informed through active acoustic telemetry, CTD casts, captive larval rearing, and field collections of related albulids and anguillids. Ocean conditions were derived from the Regional Navy Coastal Ocean Model American Seas dataset. Each spawning event was simulated 100 times using the program Ichthyop. Ten-thousand particles were released at observed and putative spawning locations and were allowed to disperse for the full 71-day pelagic larval duration for A. vulpes. Settlement densities in defined settlement zones were assessed along with interactions with oceanographic features. The prevailing Northern dispersal paradigm exhibited strong connectivity with Grand Bahama, the Berry Islands, Andros, and self-recruitment to lower and upper Abaco. Ephemeral gyres and flow direction within Northwest and Northeast Providence Channels were shown to have important roles in larval retention to the Bahamian Archipelago. Larval development environments for larvae settling upon different islands showed few differences and dispersal was closely associated with the thermocline. Settlement patterns informed the suggestion for expansion of conservation parks in Grand Bahama, Abaco, and Andros, and the creation of a parks in Eleuthera and the Berry Islands to protect fisheries. Further observation of spawning events and the creation of biophysical models will help to maximize protection for bonefish spawning locations and nursery habitat, and may help to predict year-class strength for bonefish stocks throughout the Greater Caribbean.

Introduction Connectivity of marine fish populations is driven by the complex interactions of biological processes and physical oceanographic conditions [1]. The exchange of biomass across spatially distant populations supports genetic diversity (genetic connectivity) and population growth and persistence (demographic connectivity) [2,3]. As more connections within a meta-population are formed, local populations become more resilient to perturbations and yield increased abundances [4][5][6]. Delineating population connectivity and understanding the biophysical mechanisms that drive it can improve management and conservation strategies in species for which maintenance of metapopulation connections is critical. This process is becoming increasingly important as the global fisheries paradigm moves away from discrete stock-based management towards biologically defined stocks [7]. Furthermore, delineating connectivity can allow for structured decision making when establishing spatial management strategies, such as marine reserves [8,9], challenging decisions in regions where resources are limited and enforcement is difficult [8,10].
Many marine fish species form transient spawning aggregations, migrating from their home range to disparate locations to spawn [11,12]. Identifying which local populations contribute to an aggregative spawning population is a significant component in the assessment of population connectivity [13,14]. To do so, the full extent of the catchment area-the geographic extent of adults migrating from home ranges to spawning grounds and the extent of larval dispersal [15]-needs to be established. While local population connectivity through intermixing of spawning adults is important [13,16], larval dispersal and the extent of dispersal has repeatedly been identified as the most significant driver of both genetic and demographic connectivity [1,17,18]. Dispersal is the holistic measure of connectivity, bridging spawning and recruitment events through the inclusion of transport-defined by the x, y, z displacement of larvae over time-and settlement [19]. The in situ assessment of larval dispersal source-sink dynamics in a marine environment is comparatively more difficult and labor intensive than tagging studies conducted on adult fishes [e.g., 17]. Field studies of larvae and their dispersal are often opportunistic, with region-specific goals and a broader focus on community composition [20][21][22]. Furthermore, follow up assessments of genetic relatedness [17] and/or comparing microchemical signatures with environmental concentrations [23] add time and complexity. As such, knowledge gaps often exist for information pertaining to both biological and physical controls of larval dispersal and connectivity.
Larval dispersal models (LDMs) are useful as a computational proxy for in situ monitoring of larval dispersal and movement, alleviating the logistical constraints of implementing field studies to fill knowledge gaps in larval dispersal/retention dynamics [24]. LDMs are coupled models with two modeling components, and an optional third: (1) Hydrodynamic model, (2) Lagrangian dispersal model, and (3) Individual based model (IBM). There are many regional hydrodynamic models maintained by working groups and governments, such as ROMS [25], SYMPHONIE [26], and NCOM [27]. These models numerically solve the Boussinesq hydrostatic equations of motion over a curvilinear grid to simulate the full water column ocean circulation and its hydrology. Lagrangian dispersal models build upon hydrodynamic models, implementing a particle tracking and settlement scheme for neutrally buoyant passive particles within the x, y, z, and time structure of the hydrodynamic model. Larval biology can be parameterized into the Lagrangian dispersal model, creating an IBM. IBM parameterizations include larval growth and changes in density, diel vertical migration, swimming ability, mortality, and restrictions on settlement habitat. In addition to the spatial tracking of particles over time, IBMs can also track environmental encounter histories and fate assignments [24]. Collectively, the three sub-models can be synthesized into a biophysical larval dispersal model that details the interactions between physical ocean conditions, biological processes, and settlement probability. Multiple biophysical larval dispersal modeling programs have been developed over the last two decades, such as LTRANS [28], CMS [29], Parcels [30,31], and Ichthyop [32]. The increasing accessibility and useability of these models has brought together research communities dedicated to improving simulation accuracy and ease of implementation, broadening their application, and leading to model updates as new information is discovered.
Assessing connectivity has been central to establishing and supporting the management of bonefish. Some species of bonefish, such as Albula vulpes (herein bonefish), are an economically and culturally important sportfish that supports the catch-and-release flats fisheries throughout the Caribbean Sea and northwest Atlantic Ocean [33,34]. The economic importance [35][36][37] and near threatened conservation status [38,39] of bonefish have made the species of particular interest for conservation efforts. Extensive and numerous acoustic telemetry and conventional tagging studies informed by traditional ecological knowledge have revealed much of what is known about bonefish movement, and their findings are central to how the species is managed and protected. Adult bonefish home ranges are inclusive of intertidal sand and seagrass flats, mangroves, limestone hardbottom, saltwater creeks, and channels [40]. Recreational catch-and-release fishing is focused upon the flats habitats, and release handling practices are critical to post-release survival [41][42][43]. Additionally, significant progress has been made on bonefish spawning behavior. From October through April over 4-7 day periods spanning new and full moons, bonefish migrate from their shallow-water home flats to nearshore deeper-water locations and form pre-spawning aggregations (PSAs) of 2000-5000 fish [44,45]. Currently, eight PSAs have been identified in The Bahamas [44,46] and two in Belize [47]. Home ranges that contribute to the formation of a PSA can be up to 118 km away [14,48]. From the nearshore PSA location, bonefish move offshore at dusk to spawn [44,45,49]. These tagging studies have helped to establish and support laws mandating catch-and-release fisheries [34,50], community-based management [51,52], and national park designations that protect vulnerable coastal habitats [48].
Limited information is available regarding the biology of bonefish larvae and the connectivity of bonefish populations through larval dispersal pathways. A. vulpes has an extended pelagic larval duration (PLD) and delayed settlement competency period, existing as a clear eel-like leptocephalus larvae for 41-71 days [53] before settling into littoral nursery habitats and metamorphosing into juveniles [54]. Zeng et al. [55] used an LDM for the first assessment of larval dispersal pathways for bonefish across the Caribbean, using eight known PSA locations and 18 theoretical locations. Larval dispersal restricted to the sea surface suggested significant connectivity of bonefish populations and supported the establishment of a Caribbean bonefish metapopulation for management purposes. In the time since these simulations, significant discoveries have been made regarding bonefish spawning that can better inform bonefish LDMs. Lombardo et al. [49] documented bonefish diving to depths reaching 137.9 m after moving offshore to spawn at a depth associated with the pycnocline and thermocline. Furthermore, repeated offshore tracking of bonefish from the South Abaco PSA demonstrated that the selection of a pelagic spawning site is plastic [49]. Spawning depth, the association with frontal features, and variation in selected spawning location may have a significant effect on patterns of population connectivity and year-class strength.
Here, we apply a LDM to the observed spawning events documented in Lombardo et al. [49] in South Abaco, The Bahamas for the years 2013, 2018, and 2019. The incorporation of recent discoveries in bonefish spawning behavior and the biological parameterization of leptocephalus larvae provide a more refined representation of dispersal pathways and connectivity within the meta-population. As in situ pelagic and settlement stage larval sampling within the Bahamian Archipelago has not overcome the logistical constraints of regularized sampling efforts, larval development environments are undescribed for larvae successfully dispersed to islands throughout The Bahamas. The results of the LDMs represent a baseline estimation and description of dispersal pathways, connectivity within the region, and larval development environments.
Regional settlement densities and the persistence of connectivity are used to provide suggestions for spatial management strategies directed towards preservation of population connectivity within The Bahamas. Post hoc analyses of larval development environments during dispersal allowed for the evaluation of hypothesized contrasts for larvae dispersed to different islands. While larvae inherently experience different environments throughout the dispersal process, we hypothesized that (1) central measures of environmental conditions (mean temperature and salinity over dispersal duration) and (2) dispersal pathways (dispersal linearity, settlement days post spawn, release depth, and mean dispersal depth) of successfully settled larvae maintained relative relationships between islands across years (e.g., on average, larvae that settle on Abaco have a more circuitous dispersal path than those that settle on Andros). Interannual consistency in these relative relationships would be indicative of regional prevailing conditions that normalize larval dispersal, development, and connectivity.

Spawning events
The spawning events from which the suite of larval dispersal models were created from are published in Danylchuk et al. [44] and Lombardo et al. [49], which are the only detailed documented instances of bonefish spawning movements. Offshore spawning migrations from the PSA location in South Abaco were actively tracked with acoustic telemetry in November 2013, 2018, and 2019 [See 44, 49]. In 2013, five bonefish were gastrically tagged with Vemco continuous acoustic transmitters with pressure sensors rated to a depth of 50 m (V9P-2H, 9 mm diameter, 21 mm in length, 1.6 g in air, 2000 ms transmission period). In 2018 and 2019, six and four bonefish, respectfully, were surgically implanted with Vemco continuous pressure and temperature tags with a 250 m depth limit (V9TP, 9 mm diameter, 31 mm long, 4.9 h in air, period 1000 ms; Innovasea Systems Inc., Massachusetts). Not all tagged fish were successfully tracked offshore, as some fish emigrated, or were lost due to mortality or tag failure. The numbers of fish successfully tracked in 2013, 2018, and 2019 were one, three, and one, respectfully, and in 2019 the entire aggregation was also observed by sonar during the offshore migration. Active acoustic telemetry allowed for continuous monitoring of the migration into offshore waters that exceed 1000 m, and the use of pressure sensor tags logged depth information. The spawning events were described in four-dimensions: longitude, latitude, depth, and time. In 2019, CTD casts were taken over the three days preceding the spawning event and immediately upon observation of spawning behavior using a Castaway CTD rated to 100 m (SonTek YSI, Xylem Inc., New York).
In 2013 and 2019, spawning rush behaviors-a rapid ascent rate greater than the average moving rate-and mixing movements-balling and swirling behaviors when spawning [56]were observed in the depth sensor data, indicative of the spawning events. However, gamete release over the full duration of the ascent or solely at the apex has not been confirmed. In 2018, the tag signals were lost during the initial descent in rough seas. Observations from the 2019 spawning event provided the most complete and detailed information about bonefish spawning ecology. The spawning event in 2019 revealed that bonefish spawning depth is likely associated with the pycnocline and thermocline boundary layer, as evidenced by the ejection of the acoustic telemetry tag [49]. Assuming an obligate relationship, thermocline depth was used to inform the parameterization of the 2013 and 2018 spawning events where the extent of the spawning rush below 50 m and the movement beyond the initial descent, respectively, were not observed (See Larval Dispersal Model below).

Hydrodynamic model and grid
The hydrodynamic model chosen for this study is the Navy Coastal Ocean Model American Seas (NCOM AmSeas) [27] regional model from the Naval Oceanographic Office (NAVO-CEANO). NCOM AmSeas is capable of resolving mesoscale and sub-mesoscale eddies, lending itself to applications in biophysical models [57][58][59]. The global NCOM model was developed by the Naval Research Laboratory and is based upon the Princeton Ocean Model with time invariant hybrid (sigma over z) vertical coordinates, allowing for higher resolution terrain following (sigma) and computational efficiency at deeper depths (z). NCOM AmSeas is a high-resolution derivative of global NCOM created in response to prediction and monitoring needs following the Deepwater Horizon incident [60]. The extent of the regional model covers the Gulf of Mexico and the Caribbean at 1/36˚(~3.1 km) horizontal resolution and 40 depth levels. NCOM AmSeas ingests quality-controlled observations from the Navy Coupled Ocean Data Assimilation System (NCODA) [61], consisting of satellite sea surface temperature and altimetry as well as surface and profile temperature and salinity data. Boundary conditions are applied from operational global HYCOM (1/12˚resolution) [62] and atmospheric forcings are provided by the Navy Global Environmental Model (NAVGEM) [63]. Tidal forcing is applied through two different methods. At the boundary, tides are imposed from the TPXO global inverse tidal prediction model [64]. Within the model domain, a tide-generating force is applied by incorporating astronomical forcing, ocean loading, and ocean self-attraction consistent with the TPXO model [see 65 for more detail]. NCOM AmSeas simulation archives consist of 3-hourly ocean state variables aggregated in NetCDF files including ocean temperature, salinity, sea surface height, east and northward flow velocity, and atmospheric forcing fields. Data were downloaded for each spawning season (October-April) and years (2013, 2018, 2019) observed with a spatial extent inclusive of potential dispersal to viable bonefish fisheries in The Bahamas, the northern extent of Cuba, and Florida and the Florida Keys.
All hydrography data transformations were done in MATLAB. To meet the hydrography data format and naming scheme requirements of Ichthyop, the NCOM AmSeas data were translated to a ROMS-like sigma coordinate grid. The bathymetry mask was created using the ROMS_TOOLS package [66] with bathymetry sourced from a subset of the ETOPO1 DEM [67], and clipped to only include depths greater than 1 m to eliminate land from the DEM. The shoreline mask was created from the Global Self-consistent, Hierarchical, High-resolution Geography database (GSHHG) [68]. Grid cells were masked when land comprised > 50% of the cell, single cell embayments, and single cell width edges along the domain boundaries. Spatial interpolation of the NCOM AmSeas data within the grid file was conducted using the ROMS_TOOLS package [66], which enforces mass conservation on the ROMS grid.

Larval dispersal model
We selected Icthyop [32] as our biophysical LDM program due to its frequent co-application with NCOM AmSeas data in 3-dimensional particle dispersal modeling within the Caribbean [59,69,70]. A 3 h computational timestep was used, matching the NCOM AmSeas output time interval. The computational timestep did not exceed the fine-scale limit of the Courant-Friedrichs-Lewy criterion, approximated by dt ¼ 0:7 � dGrid Umax , where dt is the timestep, dGrid is the average length of the grid cells, and Umax is the order of magnitude of the fastest current velocities in the hydrodynamic model [24,32,71]. A Forward Euler advection scheme was implemented with a horizontal dispersion rate (ε) set following equations (2 & 3) of Pelíz et al. [72]. Okubo [73] equation (4) was used to calculate the explicit Lagrangian horizontal diffusion (K h ) for a 3.1 km grid. The horizontal dispersion is randomized at each timestep along a uniform distribution (δ 2[−1,1]) [72], and the vertical dispersion is randomized using a vertical dispersion model [74] with cubic spline interpolation of the vertical diffusivity fields of the hydrodynamic model [32]. The randomized dispersion acts as a 3-D random-walk across multiple simulations, thus producing unique dispersal for each LDM iteration. To incorporate the stochastic variability of larval dispersal into the assessment of development environments, dispersal pathways, and connectivity, the LDM was run 100 times for each observed year, releasing 10000 larvae per iteration. The number of model iterations and larvae released were selected to minimize error in the potential or realized dispersal kernel, allow for representative sampling of the distributions governing the randomized directionality of the dispersion parameter enacted at each timestep, buffer against potentially significant mortality events (see Biology Parameterization), and to allow for within-year and across-year comparisons [24,75,76]. Furthermore, computational costs were considered when assessing the upper limit of how many larvae to release, as post hoc inspection would assess parameters governing every larva, at every time step, in each LDM iteration. Parameterizations for the simulation time, extent, stain (spawning event), transport, biology, and recruitment were set in the configuration XML file (see S1 File for 2019 model configuration XML).

Spawning event parameterization
The . Therefore, the observed movement and spawning behaviors were translated into proportional depths in relation to the in situ thermocline. From these proportions, the spawning movements and behaviors could be related to the NCOM AmSeas thermocline depth so that the simulated spawning event depths would be at an equally proportional depth distance from the thermocline as the observed spawning event. In 2019, the upper observed spawning depth was 63.7 m, which was 84% of the thermocline depth (75.9 m). Therefore, the upper spawning depth for each simulated spawning event would be set at 84% of the NCOM AmSeas thermocline depth. The bottom extent of the spawning depth was determined by using the observed bottom extent of the 2019 spawn and the temperature data from the expelled acoustic tag, which functioned similar to an expendable CTD capable of returning data at depths greater than the depth rating of the Castaway CTD [See 49]. The bottom extent of the 2019 spawn coincided with a deviation from the consistent temperature-depth gradient, where the rate of cooling increased with depth (i.e., the maximum first-order derivative below the thermocline). Therefore, the bottom extent of spawning for each year was placed at the depth below the thermocline where the first-order derivative was maximized in the NCOM AmSeas data. These relationships resulted in simulated spawning depth ranges of 50.36-100 m, 50.36-125 m, 41.96-150 m for 2013, 2018, and 2019, respectively. Each spawning event released 10000 particles randomly distributed throughout the water column between the bottom and peak spawning depth with a horizontal radius of 150 m; equivalent to the 300 m detection range of the acoustic telemetry tags. Larval position, salinity, temperature, and vital status were recorded every 3 h.

Biology parameterization
Little is known about bonefish leptocephalus biology and behavior throughout the pelagic duration of their development. While captively reared bonefish leptocephali displayed a

PLOS ONE
vertical orientation and sporadic short burst swimming motions [77], these observations were only during the first < 10 days of the 71-day larval period, so volitional swimming behavior throughout development is not understood well enough to inform parameterization. Although burst and critical swimming velocities of other Elopomorph leptocephali have been experimentally measured and could be used in proxy, the ecology governing sustained swimming behavior throughout ontogeny is still not known [78]. Non-volitional movement may be imparted through specific gravity, or a particle's relative density to seawater, which is likely to be important in leptocephalus dispersal as leptocephalus are uniquely buoyant [78,79]. Larval density has not been measured for bonefish to parameterize buoyancy over development. However, density measurements over development have been taken for another Elopomorph, the Japanese eel Anguilla japonica [79], which we applied to bonefish by creating an equivalency between development time and the ratio of larval density to seawater density. The measured densities of individual A. japonica were converted to a proportion in relation to the density of water at the upper thermocline. Furthermore, the length of each individual A. japonica can be converted into a proportion of the time duration of larval development. Therefore, specific gravity measures for each pre-metamorphic larva by Tsukamoto et al. [79] were converted from a length:density ratio to L=L max D l =D s , where L is the larva length, L max is the maximum observed length of the pre-metamorphic larvae, D l is the density of the larva at length L, and D s is the seawater density at the thermocline. These proportions were matched to the larval duration of bonefish (See S2 File). Diel vertical migration (DVM) in bonefish is not well documented; however, most settlement competent larvae have been caught in the top 1 m of water when light traps were set in channels exceeding 4 m deep [53]. DVM was parameterized to maintain larval position within the top 2 m of water during the settlement competency period [41-71 days post spawn (DPS)]. The DVM programming is independent of spatial proximity to settlement zones, as there is no spatial component to the DVM parameters in Ichthyop. Lethal temperatures were parameterized with a lower bound at 14˚C and upper bound at 32˚C. Congener bonefish larvae have been captured in the Gulf of California at 15˚C [80,81] and in waters up to 30˚C [21], while congener juveniles in the Indian River Lagoon, Florida have been observed in cold-snap mortality events at 14˚C [82]. Transport behavior was parameterized to a beaching particle-the computational process of particle mortality due to contacting the shore before settlement competency. The biological basis for the beaching particle was not inclusive of larvae physically contacting the shore, as coastal hydrodynamic processes would deter such occurrences. Therefore, the beaching particle parameterization would be more appropriately described as onshore mortality. Onshore mortality encompasses mortality events driven by life stage-habitat mismatches (e.g., starvation, predation, inhospitable habitat) [83,84], for which nearshore habitats inside of fringing reefs were presumed to be mismatched to pre-settlement competent bonefish larvae [53]. Additionally, the onshore mortality parameterization reduces the influence of unresolved nearshore processes over the full duration of dispersal that could not be captured at the grid resolution. Herein, the biological process of onshore mortality is synonymous with the computational process of the beaching particle parameterization.

Settlement zones parameterization
Bonefish nursery habitats are associated with littoral habitats [54]. Settlement zones (Fig 1)generalized settlement and nursery habitat equivalents that result in successful settlement if reached within the settlement competency period (41-71 DPS)-were created using ArcMap 10.6.1 (ESRI, Redlands, CA). Grid centroids were extracted from the hydrodynamic grid file in MATLAB and imported into ArcMap using the Make NetCDF Feature Layer tool. Grid points for each settlement zone were selected and converted into polygons using a 500 m buffer around the GSHHG shoreline polyline shapefile [68]. As grid points were centroids and not vertices, vertices were manually adjusted for the settlement zones to extend one to two grid cells beyond the shoreline of the land mask. Settlement zones were extended from shore to account for the potential ability of larval bonefish to maintain their position until biological and environmental conditions are suitable for settlement [85]; an ability that could not be parameterized into the biology of the larvae. Furthermore, the extension of the settlement zones ameliorates possible deficiencies in the hydrodynamic resolving of fine-scale shorelines processes. Shoreline composition was not discriminated against, as larval behavior and ability to navigate sub-optimal settlement habitats is unknown. Settlement zones from the shapefile created in ArcMap were converted to an XML file compatible with Ichthyop using MATLAB.

Assessment of dispersal
Dispersal patterns and connectivity. All statistical analyses and visualizations were constructed in R 4.0.0 (R Core Team, Vienna, Austria). Outputs from the LDMs were aggregated by year then analyzed to describe annual and prevailing dispersal patterns. Larval dispersal patterns were visualized by calculating discrete gridded larval densities at each timestep, allowing for inspection and description of dispersal pathways and the final position of each larva at the end of the competency period at 71 DPS. Fate assignments-alive, out of domain, dead beached, dead cold, and dead hot-were logged at each timestep and the mean proportion of larvae within each fate classification was calculated from the aggregated LDMs. Proportional fate assignments throughout the dispersal period were used to provide context for gridded larval densities and insight for the temporal onset of abiotic sources of mortality. Settlement clustering was quantified using kernel density estimates (KDEs) from the MASS package [86], as larvae are independent, highlighting potential hotspots for recruitment. Island connectivity was quantified for each year using the percent of LDM iterations that resulted in successful settlement, the aggregated abundances of larvae settled, and the frequency of island combinations (hereafter settlement footprints) settled by LDM iterations.
Larval development environment. The larval development environment-inclusive of temperature, salinity, depth, and measure of dispersal deflection (hereafter coefficient of dispersal)-was summarized for each successfully settled larva and assessed for variability within years in relation to which island was settled. The coefficient of dispersal was calculated as the ratio of dispersal path and straight-line distances from the spawning location to the settlement location. The closer to 0, the more circuitous the dispersal route, and the closer to 1, the more linear the dispersal route. Larval development environment characteristics of successfully settled larvae were grouped by island and year and confirmed for non-normal distributions using a Kolmogorov-Smirnov test, which indicated non-normality for all measures. A Kruskall-Wallis test was then used to evaluate whether island was a statistically significant predictor of each of the larval development environment features. Pairwise comparisons were conducted using a post hoc Dunn's test with Bonferroni correction. Pairwise comparisons were further evaluated for biologically relevant differences using Hedges' g as a standardized measure of effect size. Hedges' g is a sample size corrected measure of Cohen's d and is more appropriate when comparisons involve samples sizes � 50 [87,88]. Hedges' g is measured in relation to the pooled standard deviation of the two groups, and the strength of the effect was classified according to Cohen [89] as negligible (< 0.2), small (= 0.2), medium (= 0.5), or large (� 0.8). Confidence intervals were calculated using bootstrapping. Comparisons were conducted using the dunn.test package [90] and the effsize package [91].

Dispersal patterns and connectivity
Discrete gridded larval densities of all 100 LDM iterations for each of the observed spawning years throughout the full extent of the PLD reveal the diversity in larval dispersal pathways across years (Fig 2). In 2013, the eastern current flow through the Northeast Providence Channel forced most pre-settlement competent larvae into the east-facing southern shore of South Abaco, resulting in onshore mortality and leaving a highly reduced number of larvae to be dispersed. Remaining larvae were transported eastward into the Atlantic where most larvae were transported south. Limited larvae successfully dispersed onto Great Bahama Bank at the southern extent by transport through Exuma Sound, the Crooked Island Passage, and Mayaguana Passage. Anomalous dispersal towards Cuba was observed via the Crooked Island and Mayaguana Passages. Gyre formations north and northeast of Abaco retained relatively few larvae without indication of settlement feasibility. Larvae that left the domain did so in the southeast direction.
In 2018, larvae exhibited more diverse dispersal patterns with strong pathway consistency exhibited by high gridded densities. Again, pre-settlement competent larvae experienced onshore mortality along the east-facing southern shore of South Abaco, though magnitudes less than in 2013, but also along the south-facing shore of South Abaco, resulting in onshore mortality. Larvae that did not experience onshore mortality proximal to the spawning location were mostly retained by eddies within Northeast Providence Channel and the large ephemeral gyre within Northwest Providence Channel. Strong northern transport from both Providence Channels was observed, with a substantial abundance of larvae exiting the model domain to the north-northwest, most of which originated from Northeast Providence Channel. Increased abundances on the north side of Grand Bahama indicate the presence of a larvae retaining eddy feature. Lesser southern transport was observed into the Tongue of the Ocean (TOTO), and fewer into the Exuma Sound. Eddy features retained larvae at both the northern and southern extents of the dispersal region, though it did not appear that larvae returned to the Bahamian Archipelago once circulated beyond its extent. Westward dispersal extended out to Bimini and was limited by the Gulf Stream, while eastward dispersal was controlled by Antilles Current dynamics. Few larvae were dispersed to the northern shore of Cuba via TOTO and traversing over the shallows of Great Bahama Bank.
The 2019 dispersal patterns were most similar to the patterns seen in 2018, with strong retention within the Providence Channels and large numbers of pre-settlement competent larvae experiencing onshore mortality along the south-facing shore in South Abaco. The majority dispersal pattern interacted with the large ephemeral gyre positioned within Northwest Providence Channel, with larvae transported westward out of Northwest Providence Channel and north along the Gulf Stream, exiting the model domain to the north. Larvae exited the domain over a vast area bounded to the west by the Gulf Stream and uninhibited to the east. Larvae did disperse south into TOTO and south along the Antilles Current, though their southern extent was limited to the eastern nearshore region of Andros and to the eastern nearshore region of Long Island. The western extent of dispersal reached Bimini and continued further west beyond its shore, and larvae uniquely occupied the area between Bimini and Andros. Similar to 2018, abundances north of Grand Bahama indicate an established eddy. In contrast to 2018, transport northward out of Northeast Providence Channel was not as significant as that from Northwest Providence Channel. Larval transport out of the domain was also more diffuse, and a gyre northeast of Eleuthera reduced the strength of transport south along the eastern shore of Eleuthera into Exuma Sound.
The final position of all larvae across aggregated LDM iterations in 2013 shows the highest density of larvae were found along the southern shore of Abaco, with 95% of larvae within one grid cell (Fig 3) as a product of pre-settlement competent onshore mortality (Fig 4). Larvae dispersed eastward beyond the model domain's 72.5˚W limit, with most leaving the domain southeastward. No apparent retention features were visible given the absence of clustering. In 2018, 46% of larvae were found within one grid cell along the southern shore of Abaco as a result of pre-settlement competent onshore mortality. Many larvae remained within the Northeast Providence Channel after 71 d, with larvae retained by the eddies within Northeast Providence Channel and the ephemeral gyre within Northwest Providence Channel. Larvae continued to be projected from the Northeast Providence Channel eddies in a northward trajectory, and a cluster of larvae were positioned within the shallow waters of Great Bahama Bank north of Andros. The majority of larvae left the model domain through the northern boundary. In 2019, 14% of larvae were found within one grid cell along the south-facing shore of South Abaco, mostly as a result of pre-settlement competent onshore mortality. Similar to 2018, larvae were retained within the Northwest Providence Channel ephemeral gyre, though in much greater numbers. Larvae continued to be projected from the gyre northward along the Gulf Stream. Larvae were also similarly clustered north of Andros, but over a much larger area spanning 100 km across. Large numbers of larvae were diffusely distributed throughout the sea northeast of the Bahamian Archipelago, and larvae exited the model domain over the full expanse to the north and northeast. The final positions of larvae across years were mostly determined by pre-settlement competent onshore mortality events (Fig 4). All years experienced � 50% pre-settlement competent onshore mortality within the first day, with 95% pre-settlement competent onshore mortality in 2013. No heat mortalities were experienced in any year and cold mortalities were observed in 2018 and 2019, accounting for 4.7% and 4.6% Larval settlement was successful for < 0.01% of larvae in the aggregated 2013 LDM iterations, with five island groups settled upon. Settlement KDEs indicate that the core settlement distribution was located at the southern extent of The Bahamas, at Acklins/Semana Cay ( Fig  5). In the aggregated 2018 LDM iterations, 0.12% of larvae successfully settled upon nine island groups with core settlement distributions located at the north end of Andros and west end of Grand Bahama. In the aggregated 2019 LDM iterations, 0.38% of larvae successfully settled upon eight island groups with core settlement distributions located at the southern end of Abaco, the Berry Islands, and mid-longitude of Grand Bahama.
Connectivity between Abaco and other islands varied across years, with the strongest connectivity showing island-level self-recruitment to Abaco and strong larval supply to Grand Bahama (Fig 6). The relationship between the number of larvae settled and the number of LDM iterations settled per island was positively monotonic. Connectivity in 2013 was especially poor given the high onshore mortality. The percent of 2013 LDM iterations which yielded successfully settled larvae was low, with larvae settling upon Acklins/Semana Cay most frequently at 14% of LDM iterations, though larval abundance was low. The most common 2013 LDM iteration settlement footprint-combinations of islands settled-occurred in 14% of the iterations and consisted of only Acklins/Semana Cay (S2A Fig). In 2013, no LDM iteration yielded more than 0.02% successfully settled larvae, which occurred in 4% of the LDM iterations, and the median percent of successfully settled larvae per LDM iteration was 0%.
Connectivity in 2018 was substantially stronger, with more consistent and abundant larval sourcing. The percent of LDM iterations settled in 2018 was nearly 100%, with Abaco settled most frequently at 96% of LDM iterations. Successful settlement in 2018 decreased with latitudinal distance from Abaco, as did abundances, though Eleuthera and New Providence had much lower settlement success than Andros. However, this could be attributed to differences in shoreline extents. Islands south of New Providence showed connectivity of little importance. The most common settlement footprint, at 22% occurrence, was Abaco, Andros, and the Berry

PLOS ONE
Islands. In 2018, the median percent of successfully settled larvae was 0.1%, which occurred in 5% of LDM iterations, while the minimum of 0.03% occurred in 3% of iterations and the maximum of 0.23% occurred in 1% of iterations.
The 2019 spawn supplied the highest abundances and consistent sourcing of larvae. Four islands-Abaco, Grand Bahama, the Berry Islands, and Eleuthera-had larvae settle in 100% of LDM iterations, with both Andros and New Providence above 50%. The most common settlement footprint occurred in 32% of LDM iterations, included Abaco, Grand Bahama, the Berry Islands, Eleuthera, Andros, and New Providence. In 2019, the median abundance of successfully settled larvae was 0.37% and occurred in 3% of LDM iterations, while the minimum of 0.17% successfully settled larvae occurred in 1% of iterations, and the maximum of 0.58% successfully settled larvae occurred in 2% of iterations.

Larval development environment
Larvae were grouped by year and island settled to quantify the variability in development environments experienced and to assess whether relative relationships are consistent across years.

PLOS ONE
Due to the nearly complete onshore mortality event of larvae from the 2013 spawn, statistical comparisons were only possible between the 2018 and 2019 spawns. Across years, larvae that settled more proximal to the spawning location had lower coefficients of dispersal, indicating retention in circulation features for longer distances (Fig 7). Successfully settled larvae ranged in their coefficient of dispersal from < 0.0001 to 0.620, equivalent to dispersal paths that were >1000 times to approximately 1.5 times the straight-line distance from the spawning location

PLOS ONE
to the settlement location. The median coefficients of dispersal across years were 0.47, 0.13, and 0.07, respectively, with median dispersal path distances of 986.1 km, 1050.5 km, and 1057.5 km. Island specific coefficient of dispersal distributions were similar in shape and range in 2018 and 2019, which is indicative of similar or prevailing hydrodynamic conditions. However, there were differences in the coefficient of dispersal distributions for Abaco and Grand Bahama between years, as the gyre and eddy structures near the spawning locations were different in position, size, and numbers (S2E and S2F Fig). Larvae exhibited within year island level differences in coefficients of dispersal (S1A and S1B Table), with effect size differences in 2018 smaller than those in 2019 (Fig 8). Of the total 43 contrasts made within 2018 and within 2019, 36 had large differences with the difference in means � 0.8 of the pooled standard deviation (sd) per comparison. Relative relationships between islands were largely conserved between the 2018 and 2019 spawns, with 73% of relationships maintained, indicating stable dispersal pathways.
Settlement most frequently occurred within the first third of the settlement competency period for all years, with settlement occurring throughout the 41-71 d settlement competency period (Fig 7). reached the 50 th percentile of their total LDM aggregated abundance later than those above this latitude. Larvae exhibited within year island level differences in settlement DPS (S1C and S1D Table), with effect size differences in 2018 smaller than those in 2019 (Fig 8). Large differences in effect size were only seen within the 2019 spawn, with 10 comparisons having a difference in means � 0.8 of the pooled sd per comparison. The large differences in 2019 were mostly comprised of comparisons with Bimini, where larvae dispersed to Bimini settled earlier than other islands. Relative relationships between islands were mostly conserved between the 2018 and 2019 spawns, with 53% of relationships maintained.
Larvae were released randomly throughout the parameterized spawning depths for each year and successfully settled larvae were produced from nearly the full range of depths. Across all years, larvae were more likely to successfully settle if released at a deeper depth (Fig 7). However, successful settlement did diminish at the bottommost extent of the spawning depths, which was parameterized to coincide with the depth below the thermocline where the temperature-depth gradient grew. Very few larvae released above the thermocline successfully settled in any year. Island was not significantly associated with release depths within either year (S1E and S1F Table), and no comparative effect sizes exceeded a Hedges' g of 0.8 (Fig 8). The mean release depth of larvae dispersed to each island within years was statistically equivalent, with only seven of the 43 contrasts having confidence intervals that did not overlap with zero, most of which were comparisons with larvae that settled upon Bimini.
Larvae successfully dispersed over a wide range of depths across years, apart from the few larvae that settled from the 2013 spawn (Fig 7). Most successfully settled larvae were dispersed at a mean depth at or near the thermocline, though having been released below the thermocline and often much deeper. The distribution of mean dispersal depths in 2018 was broader than the distribution for the 2019 spawn, in part due to the heavily skewed shallow mean dispersal depths for larvae settling upon Andros in 2018. While mean dispersal depths were statistically different in within year, across islands comparisons (S1G and S1H Table), the mean dispersal depths only differed in effect size in 2018 when comparisons were made with Andros (Fig 8). As such, mean dispersal depths exhibited stable relative relationships, in that mean dispersal depths did not differ.
Mean salinities experienced by larvae over their dispersal did not vary more than approximately 0.6 PSU. Minimum and maximum salinities were also examined, though they collectively varied by only 2 PSU. Such minute differences given the euryhaline classification of Albulid larvae [92] did not warrant further investigation of differences in salinity exposure during the dispersal process or potential implications. Mean temperatures experienced by successfully settled larvae ranged from 21.73˚C to 26.40˚C (Fig 7), the coldest minimum temperature approached within 0.02˚C of the 14˚C cold mortality limit, and the warmest temperature was 28.95˚C. Across all years and all larvae, the average range in temperatures experienced was 9˚C, with a minimum and maximum deviation of 3.5˚C and 14.4˚C, respectively. The coldest mean temperatures were experienced by larvae dispersed to Abaco and Grand Bahama, with larvae experiencing warmer waters when dispersed to more southern locations. Mean temperatures were similarly ranged in 2018 and 2019, with the medians approximately at 25˚C. Within year island comparisons were statistically different (S1I and S1J Table), and 15 of 43 contrasts having large differences in effect size (Fig 8). Differences were larger in 2018, though 67% of relative relationships maintained their directionality despite smaller differences among islands.  selection via these cues have been observed in other pelagic spawning fishes, such as bluefin tuna Thunnus thynnus [93], blue whiting Micromesistius poutassou [94], and many eel species that share an evolutionary history with bonefish [95]. Further efforts in bonefish active tracking and hydrographic sampling are needed to confirm these hypotheses.

Dispersal patterns and connectivity
Mortality. Flow direction, speed, and proximity to shore at the time of spawning significantly influenced the proportion of larvae available to settle. Pre-settlement competent onshore mortality within hours of the spawn was observed across all years, driven by spawning in close proximity to the western shore of South Abaco in the presence of an east flowing current. The significant and nearly instantaneous mortality may have been reduced given a finer grid system that could capture nearshore hydrodynamics; however, this LDM framework could not accommodate the nested grid systems that would be required, nor were such efforts feasible given computational costs and the size and complexity of the model domain. Furthermore, the 1/36˚grid resolution already exceeded the resolution-domain relationships of most biophysical models within the literature [24]. The cessation of pre-settlement competent larval dispersal upon programmatic beaching/onshore mortality inherently reduces the influence of unresolved nearshore processes-inclusive of abiotic (hydrodynamics) and biotic (ecology)by preventing continued or repeated exposure to nearshore regions of uncertainty. In doing so, significant losses of larvae may and did occur-the 2013 LDM resulted in 95% onshore mortality across all iterations. It is important to highlight the limitation of the LDM construct in that it is not feasible to parse the abiotic, biotic, and programmatic fractions of onshore mortality. Thusly, a tradeoff of uncertainty in dispersal for uncertainty in the source and true extent of onshore mortality is chosen when onshore mortality in enabled. Estimations in dispersal and environmental exposures will better align with the resolution and processes in which the hydrodynamic model was structured to capture. However, as this is applied to bonefish larvae, the absence of pre-settlement competent larvae in nearshore larval sampling where Albula spp. spawn supports that the pre-settlement competent larvae would not survive within region where nearshore hydrodynamic and ecological processes exert control on dispersal viability [20,21,53,80], and thus the programmatic application of onshore mortality is suitably applied.
While these LDMs programmatically implicated physical contact with the shore as a cause of mortality to larvae, this is unlikely and without documentation. However, the mechanisms for mortality upon encountering a habitat before settlement competency are well documented. Habitat mismatch with the onset of endogenous feeding or thermal thresholds is a significant source of larval mortality [84], as is predation by planktivores in both pelagic and reef environments [83] that fringe the shores throughout The Bahamas. Furthermore, nearshore turbulence would likely lead to mortality in underdeveloped, fragile leptocephalus. Beyond thermal and domain boundary mortality, the multitude of sources for mortality in early-stage larvae were not possible to account for and parameterize within this LDM framework. Pre-settlement competent onshore mortality did appear to result in mortality at rates similar to other models that did incorporate mortality caused by mismatches in prey availability and predation [96], therefore the programing parameterization of beaching/onshore mortality may be useful as a proxy for mismatch driven mortality when parameterizations are unknown or not readily implemented into the LDM framework. Further work on the feeding ecology of bonefish larvae needs to be done to identify prey sources throughout the entire larval stage, as this information is unknown and would be valuable to incorporate into future LDMs as a spatiotemporal dispersal viability restriction.
The thermal tolerances of bonefish larvae are not well documented and can only be implied by capture data. We estimated the lower limit of thermal tolerance to be 14˚C as informed by capture data of Albula sp. within the Gulf of California [80,81] and mortality of juveniles within the Indian River Lagoon, FL, USA [82]-neither of which are A. vulpes. Cold mortality events within the 2018 and 2019 LDMs occurred early during dispersal at the initial interaction with the ephemeral gyre within Northwest Providence Channel. The cyclonic gyre subducted larvae from the deepest extent of the spawning range into deeper waters at approximately 250 m that reach 14˚C. Bonefish larval dispersal and early life history is tightly coupled to interactions with gyre and eddy features as evidenced by their extended pelagic larval duration. The observation of gyre mediated thermal mortality for larvae spawned at the deepest extent of spawning rushes introduces a mechanism for limiting the depth at which spawning rushes begins for bonefish. Behavior implicated in sensing the boundary of preferred spawning habitat was observed through acoustic telemetry data, where bonefish oscillated across depths, seemingly to inspect hydrologic conditions throughout the offshore migration process [49]. Such behaviors allow bonefish to target the thermocline at the upper limit of the spawning rush and also appears to appropriately limit the bottom extent of the spawning rush to reduce mortality, as evidenced by the reduction in larval settlement success at the bottom most extent of the spawning rush (Fig 7; S2G Fig). Conversely, no upper thermal mortalities occurred in any year. However, the upper thermal tolerance for bonefish larvae should be investigated given predicted climate scenarios that will increase heat retention in surface waters and potentially destabilize current patterns and winter stratification [97].
Dispersal pathways. Larvae that did not experience mortality shortly after hatching circulated up to 71 DPS, allowing for inspection of unique and prevailing hydrodynamic conditions, including those associated with increased settlement success. Dispersal patterns across the years can be classified into two paradigms: Northern dominant and Southern dominant (Fig 9). Northern dominant dispersal is supported as the prevailing dispersal pattern given the consistency between the 2018 and 2019 LDMs. Furthermore, the underlying hydrodynamics follow known prevailing seasonal circulation patterns [98], and regional LDMs from other studies and species show dispersal patterns consistent with the Northern dominant paradigm [55,99].
We did observe variability within the Northern dominant dispersal paradigm. In 2018, the presence of clustering mesoscale eddies within Northeast Providence Channel forced diffuse dispersal as portions of the larvae were divided among the periphery of a cyclonic eddy against the west facing shore of South Abaco and two anticyclonic eddies to the south and east of South Abaco. The anticyclonic eddy to the east subsequently passed a majority of the larvae to the north-flowing Antilles Current that flows along the eastern shores of Abaco, which was the main source of larval dispersal out of the Bahamian Archipelago and the model domain. Larvae within the southern anticyclonic eddy were moved south into TOTO and dispersed across New Providence and eastern Andros. The cyclonic eddy pushed larvae into the Northwest Providence Channel ephemeral gyre from the south, as the eddy and south flowing currents prevented larvae from dispersing westward immediately post spawn like what was observed in 2019. The 2019 spawn yielded the most successfully settled larvae and had the smaller of the two settlement footprints. A majority of the larvae were immediately drawn into the Northwest Providence Channel ephemeral gyre, mediated by the presence of an anticyclonic eddy positioned more proximally to the Little Bahama Bank between Abaco and Grand Bahama. Undisrupted inclusion into the Northwest Providence Channel gyre generated a slow and consistent dispersal of larvae to islands that interface with both Providence Channels. This also passed larvae westward into the Gulf Stream where larvae were lost from the northern domain.
The presence of the ephemeral gyre within Northwest Providence Channel undoubtedly is the most important circulation feature in determining the larval dispersal dynamics of bonefish sourced from the South Abaco spawn. Its presence provides a prevailing NW dispersal distribution and promotes larval retention within the northern extent of the Bahamian Archipelago. The regularity in which the gyre forms over winter [98] provides another predictable hydrodynamic feature in addition to the predictably stable thermocline that bonefish can rely upon to maximize larval survival [100]. Predictable and prevailing environments are central to the migratory spawning strategy, as significant energy investments are made prior to undergoing migrations to distant locations where conditions are unknown [11,101].
The Southern dominant dispersal paradigm appears to be anomalous, only strictly occurring in 2013. Successful settlement was low given this paradigm, providing more support for the Northern dominant dispersal paradigm as the prevailing condition. However, the difference among the two paradigms might not be as large as presented here if not for the near total onshore mortality event in 2013. Conversely, there is the potential for complete mortality in

PLOS ONE
the absence of prey availability if sufficient prey densities for larvae are reliant upon the presence of aggregating features, such as gyres and eddies [102], that were notably absent in the 2013 dispersal paths. The Southern dominant dispersal paradigm is characterized by an interaction with the Antilles Current, in particular the deep waters and counter current. The surface currents of the Antilles Current flow in a northward direction, while the counter current and deep waters flow south [98]. Evidence of a sustained southward flow along the eastern boundary of the Bahamian Archipelago can been seen across all years, with larval deposition varying between Eleuthera and Acklins. Strict adherence to the Southern dominant dispersal paradigm is anomalous, but the underlying currents do appear to be regularly involved in shaping the full dispersal footprint of the Northern dominant dispersal paradigm as well given the presence of southern dispersal along the eastern margin of the Bahamian Archipelago in 2018 and 2019.
These LDMs were uniquely constructed to incorporate new empirical information to simulate outcomes from observed spawning events. However, this effort was not the first to examine bonefish larval dispersal through LDMs, as Zeng et al.
[55] did so to estimate regional connectivity of the Caribbean bonefish metapopulation. Comparisons between the two model structures and time periods allows for a more longitudinal assessment of the dispersal pathways, though parameterizations and metrics of success do vary considerably (Table 1). Zeng et al.
[55] summarized dispersal as the final position for larvae at day 53, which is comparable to this study (Figs 2 and 3), as their final positions reported are not truly settlement. The density of final positions appears to capture the overall variability in dispersal observed in the LDMs within this study, with dispersal throughout the entire Bahamian Archipelago and Cuba; however, vagrants from the Bahamian Archipelago to the north are fewer in relative abundance to the rest of the domain and dispersal to Cuba is substantially greater. These differences can be attributed to contrasting programming in addition to hydrodynamic variation, inclusive of yearly, monthly, dispersal depth, and resolution differences. The density of final and this study that the prevailing dispersal pattern is that of the Northern dominant paradigm and that larvae are most consistently distributed-independent of settlement-within Northeast Providence Channel. Connectivity. The two larval dispersal paradigms yielded substantially different connectivity networks, with the Southern dominant dispersal paradigm in 2013 uniquely sourcing larvae southeast of San Salvador and beyond the southeastern extent of The Bahamas. In both paradigms, spawning events in Abaco have poor connectivity to islands southeast of New Providence. The Northern dominant paradigm sustained strong levels of self-recruitment and connectivity with islands north of New Providence, including Andros. Given that the hydrodynamic conditions in 2018 and 2019 align with prevailing conditions, the resultant settlement footprints can be inferred as the prevailing connectivity network for bonefish larvae sourced from Abaco. The prevailing connections are to Abaco, Grand Bahama, the Berry Islands, Eleuthera, Andros, and New Providence.
The viability of these connections is dependent upon the environments in which larvae settle, as habitats are not always suitable for a given life stage, disrupting the completion of connectivity (i.e., next generation spawning) [103]. We did not discriminate against shoreline composition, which can preclude larval settlement, especially for bonefish that require softly sloping sandy shorelines with low energy littoral zones [104]. Settlement positions within the LDMs are not necessarily definitive given that bonefish, like other fish [85,105,106], are likely to exhibit some capacity of spatiotemporal control during settlement to select for compatible nursery habitat. At Lee Stocking Island, bonefish recruitment showed strong periodicity with lunar or tidal cycles [53], suggesting that larvae can maintain their position until nighttime flood tide conditions are conducive to safe passage over the reefs [83] and into nearshore nursery habitats. The duration and conditions in which bonefish can maintain and delay settlement is unknown, and not a behavior that was possible to parameterize within the Ichthyop LDMs. We made an effort to capture the ability of positional maintenance within the settlement competency period by extending the settlement zones one to two grid cells beyond the shoreline (7.4 km maximum distance from shore); however, doing so may not entirely capture the behavior if larvae are able to maintain position before the 41 d lower limit of the settlement competency period. Given the unknown positional maintenance behaviors and the programmatic limitations in accounting for such behaviors, the assumption that successful settlement in the LDMs is equivalent to successful recruitment is required.
Post hoc ground truthing in future studies would allow for scenario-specific validation of LDM parameterization and results. Concurrent studies of genetic relatedness have proven useful in confirming connectivity and suggesting management strategies [70]. Assessments of genetic relatedness and connectivity independent of LDMs have also been applied to species in The Bahamas, including bonefish [107]. Douglas et al. [107] analyzed bonefish relatedness between Eleuthera and Grand Bahama, finding that bonefish on Grand Bahama's ocean side (i.e., Northwest Providence Channel) and north exposed flats side are sink populations, as are bonefish at the southern extent of Eleuthera. These results and the results presented in seining efforts throughout The Bahamas can be used to provide context and ground truthing to the prevailing conditions observed in the LDMs within this study. Both Grand Bahama locations were strong recipients of larvae from the 2018 and 2019 spawns, while no larvae settled in the Eleuthera locations identified as sinks. While Douglas et al. [107] observed south-north source-sink dynamics between the Eleuthera and Grand Bahama, it is within reason that prevailing conditions promote similar dynamics between Abaco and Grand Bahama in support of the LDM observations. On Eleuthera, LDM larvae settled along the northern and eastern shoreline, while few larvae settle upon the leeward shores. Rigorous sampling for juvenile bonefish has only been conducted on Eleuthera, and efforts support prevailing transport to the windward side of Eleuthera where protected embayments contain nursery habitat [104]; however, CPUE of juvenile bonefish were comparatively higher at leeward stations within Rock Sound. Larvae in the 2018 suite of LDMs did reach Rock Sound in higher relative abundances than the windward embayments, though arrival to settlement habitat preceded the lower end of the settlement competency period (see Fig 2). Additional simulated sourcing of larvae to locations where juveniles have been found include Fresh Creek of Andros [108], Fish Creek and the Marls of Abaco (J. Lewis, pers. comm.), and August Creek and Pelican Creek of Grand Bahama (J. Lewis, pers. comm.).

Larval development environment
Environments exert influence on physical performance characteristics during larval dispersal and development, as larvae exhibit phenotypic plasticity in traits such as growth rate, thermal tolerance, osmoregulation, and swimming performance [109]. Variability in phenotypic expression does influence the true measure of connectivity as phenotype-environment mismatches can reduce individual fitness, preventing the full cycle of connectivity where an individual further contributes offspring to the population [110]. Accordingly, assessing the larval development environments experienced by larvae as they disperse throughout The Bahamian Archipelago was important for gaining insight into significant variability that could confer advantages or deficiencies to larvae dispersed to certain islands.
The relative relationships of larval development environments for successfully settled larvae were largely consistent between 2018 and 2019, as both years followed the Northern dominant dispersal paradigm. However, the magnitude of differences was not equal between years, which can be attributed to variation in the hydrodynamics. Conversely, many comparisons indicated no differences between islands. Collectively, consistent relative relationships and null comparisons are indicative of prevailing conditions and consistency in the larval development environments across space and time.
The strongest contrasts in larval development environments per island were in the coefficients of dispersal and settlement DPS (Fig 7). The two metrics were not correlated, though both could potentially have strong influences on larval condition and success in adapting to littoral nursery habitats. Many of the islands in which larvae successfully settled had small coefficients of dispersal, implying that their dispersal pathways were heavily influenced by interacting with circulation features. Beyond their larval retention capacities, gyres and eddies create environments that have unique physical properties that strongly influence both bioticabiotic and biotic-biotic interactions [111].
The primary circulation feature in the northern extent of the Bahamian Archipelago is the mesoscale gyre found within Northwest Providence Channel. The gyre forms in winter and is created by an influx of cool Atlantic water from the east through Northeast Providence Channel and warm Gulf Stream water from the west through Northwest Providence Channel [98]. Inputs from the Atlantic to the east and the Gulf Stream to the west are not equal or always consistent in their volume transport [98]. In 2013, this resulted in a deeper mixed layer depth. In 2018, a warm core anticyclonic eddy was present due to elevated velocity intrusion of cold water from Little Bahama Bank [98]. And in 2019, a slower velocity warm core anticyclonic eddy was present (S2E- S2G Fig). Gyre presence under prevailing conditions places larvae into highly productive conditions, where both upwelling from cyclonic eddies within the > 1000 m deep Providence Channels and the downwelling anticyclonic gyre can supply nutrients for phytoplankton and zooplankton to thrive in an otherwise nutrient deficient environment [112,113].
The stable thermal stratification of the region in winter and eddy presence appears to be a larval development environment that is tightly coupled to bonefish evolutionary history as congeners in the Gulf of California have been found in their highest abundances within the same environment [21]. Accordingly, larvae with higher coefficients of dispersal may have increased access to prey, and subsequently have improved body condition and higher survival upon settlement. The positive relationship between body condition and settlement survival has been documented in reef fish [114], and body condition can even offset deficiencies in competitive ability brought on by late settlement (i.e., priority effects) [115][116][117]. Bonefish larvae settling in The Bahamas are likely to be less dependent upon the compensatory competitive advantage of settling with higher body condition than reef fish due to nursery habitat composition and obligations. Bonefish exploit unstructured littoral nursery habitats that are not constrained by competitive or predatory priority effects, as such habitats in The Bahamas are plentiful and the success of their recruitment appears to be facilitated by inter-species schooling with gregarious small-bodied fish [118,119]. More research is needed to understand how body condition in settlement stage bonefish effects recruitment.
Differences between the other successfully settled larval development environment measures-depth released, mean dispersal depth, and mean temperature-were much smaller and more inconsistent in their relative relationships (Figs 7 and 8). However, consistency across islands can support the presence of predictable and stable environments that are vital to larval survival. The peak distributions of mean dispersal depths were shallower than release depths across all islands and were closely associated with the thermocline. Larvae are universally found more frequently associated with stratified oceanographic features and have both biological and behavioral means to maintain close association [112]. The boundary layer within the stratification feature provides insular protection from the larger stochastic processes that occur within the larger water masses, providing reduced turbulence and flow velocities [120] over large amplitude orbital motions [100], more predictable dispersal [121], and prey availability [112]. The strong association of the deep chlorophyll maximum and the thermocline within the eastern Caribbean exemplifies the nutritional benefits of thermocline associations for larval fish within The Bahamas [122]. The buoyancy of bonefish eggs and larvae [123,124] are most likely to have imparted the consistent relationship between the thermocline and mean dispersal depth. An instantaneous buoyancy rate of 8.5 m�d -1 at 0 h post-hatch (25 h post-fertilization) was calculated using the standard Stokes' Law following the methods in McDonnell and Buesseler [125] with a larval density of 1.02 × 10 −5 kg�m -3 , head diameter of 2 × 10 −4 m, a length of 0.005 m [123], a dynamic viscosity of 0.959 kg�m -1�s [126], and a water density of 1024.4 kg�m -3 . The buoyancy rate allowed for larvae to rise towards the thermocline, and the rate of decreasing buoyancy over development reached a nearly neutral specific gravity in close temporal proximity to the larvae reaching the thermocline.
Temperature also significantly impacts larval development, growth, and survival [109]. While temperature effects on bonefish larvae are unknown, Japanese eels in captivity experienced higher hatching and survival rates, as well as faster hatching, prefeeding growth, and onset of feeding in waters progressively warmer from 19-28˚C [127]. Similar results were observed in tropical reef fish [128]. Consistency in the positive relationship between temperature and development suggests warmer conditions would benefit bonefish larval development; however, the lower and upper range of beneficial temperatures cannot be assumed. Future research into how temperature impacts bonefish larval development is needed in order to fully assess the implications of variable temperatures during dispersal.

Future management and bonefish research suggestions
Zeng et al. [55] concluded that management plans effectuated at scales pertaining to isolated stocks were not sufficient to effectively conserve bonefish populations within the Caribbean. The same issue of scaling can be applied at an intra-jurisdictional scale within The Bahamas. While the harvest of bonefish is restricted to subsistence fishing within The Bahamas [33], loss of PSA and nursery habitats remains a present threat to bonefish populations [38]. The Bahamas has established 77 marine protected areas (MPAs), marine managed areas (MMAs), fisheries reserves, ecological reserves, and national parks that protect marine and shore habitats [129], with bonefish as a focal species for habitat protection in some of these areas. Tagging studies focusing on the movements of bonefish from home flats to PSA locations [14,48] helped to inform park inclusion of bonefish habitats. Data on nursery habitats and settlement patterns did not exist for consultation and incorporation into spatial management plans. The results of our study show an alignment of core settlement areas and park protected areas throughout Abaco, Grand Bahama, the Berry Islands, and Northern Andros. Within these locations, protections from development and water quality degradation should be sufficient to protect nursery habitats and larvae sourced from spawning events in Abaco. Our models did not discern potential bonefish nursery habitats through remote sensing or in situ assessments. All shorelines were included as potential settlement zones, as the capabilities of bonefish larvae to move through habitats is unknown. Furthermore, methods to delineate bonefish nursery habitats are still being refined and are labor intensive [54,104]. Accordingly, suggestions to include other locations with high concentrations of larval influx to the Bahamas' parks and spatial protections programs requires assessments of habitats to delineate nursery status. We suggest further examination and consideration of potential nursery habitats in the core settlement areas of Moore's Island Abaco (northwest of the spawning locations), the central most eastern islands of the Berry Islands, and North Eleuthera that would result in the creation of new parks. Expanding and conjoining of Marls of Abaco National Park (west central Abaco) to Cross Harbor National Park (South Abaco) would encompass the 23 km of catchment to the north that contributes to the South Abaco PSA [14], and thus would be inclusive of self-recruiting larvae. Larval supply was strong along the northern shore of Grand Bahama, and westward expansion of North Shore/The Gap National Park (north central Grand Bahama) to conjoin with West End MPA (northwest Grand Bahama) would protect recruitment to this region. A core settlement area was observed upon the northern tip of Andros, where Joulter Cays National Park encompasses the sand flats habitats that North Andros bonefish occupy. However, Joulter Cays National Park does not extend to the creeks and marls between the sand flats and mainland of North Andros -habitats that are likely to support juvenile bonefish recruitment. Expansion of Joulter Cays National Park southward to include the creek and marl habitats is recommended. Core settlement areas were also observed in East New Providence and West Grand Bahama; however, both regions are heavily developed by nearby Nassau and Freeport, which have already contributed to the destruction of nursery habitats and low bonefish densities found in those regions.
As new biological information is discovered on the spawning and early life history of bonefish, the construct of bonefish LDMs should be further refined to increase accuracy in representing the dispersal of larvae. We suggest future research subjects that will contribute to bonefish LDMs. First, prioritize describing ontogenetic changes in larval swimming capabilities, in both the pelagic and nearshore environments, since larval behavior can significantly influence dispersal and connectivity [76]. Second, evaluate the relationship between bonefish spawning site selection and physical ocean conditions. Additional bonefish spawning tracks should be conducted while implementing the use of a CTD with a depth rating > 100 m to capture the bottom extent of the offshore migration in instances where fish move below 100 m. The required equipment will be a logistical constraint, as the CTD and retrieval equipment needed to sample such depths increase in size and mass. Finally, future bonefish LDM studies should incorporate ground truthing using juvenile habitat surveys [see 105] and/or larval collections [see 53] and evaluate genetic relatedness. Collectively, these research objectives will fill knowledge gaps in bonefish early life history that have persisted for decades, and ultimately help to inform conservation measures across the Greater Caribbean.